function [gbest,gbestval,fitcount] = FCO_func(fhd,Dimension,Particle_Number,Max_Gen,lowerbound,upperbound,pos,bool_Plot,varargin)

gbest = 0;
gbestval = 0;
fitcount = 0;
rand('state',sum(100*clock));
me = Max_Gen;
ps = Particle_Number;
D = Dimension;
g = 0.13;

num_FClevel = 6;
prob_forage = 0.9;
prob_invdir = 0.5;
dir_forage = 1;
size_each_level = repmat(floor(ps/num_FClevel),1,num_FClevel);
size_each_level(1) = ps - (size_each_level(1)*(num_FClevel-1));
% for n = 1 : ps
%     if(n<=ps-size_each_level(num_FClevel))
%         if(mod(n,size_each_level(1)) == 0)
%             level_species(n,1) = n/size_each_level(1);
%         else
%             level_species(n,1) = floor(n/size_each_level(1))+1;
%         end
%     else
%         level_species(n,1) = num_FClevel;
%     end
% end
lb_level(1,1) = 1;
ub_level(1,1) = size_each_level(1);
for l = 2 : num_FClevel
    lb_level(l,1) = sum(size_each_level(1,1:l-1))+1;
    ub_level(l,1) = sum(size_each_level(1,1:l));
end

[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);
% pos = VRmin+(VRmax-VRmin).*rand(ps,D);

fitcount = ps;
vel = zeros(ps,D); % initialize the velocity of the particles

fitness = feval(fhd,pos',varargin{:});
[fitbest,gbestindex] = min(fitness);
fitworst = max(fitness);
gbest = pos(gbestindex,:);
gbestval = fitbest;
if(bool_Plot)
    figure(8)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('FCO');
end
% return
for iter = 2:me    
    if(fitbest == fitworst)
        mass = repmat(1/ps,ps,1);
    else
        mass = (fitness'-fitworst)/(fitbest-fitworst);
    end
%     mass = mass / sum(mass);
    sort_mass = [mass (1:ps)'];
    sort_mass = sortrows(sort_mass,1);
    
    a=2-iter*((2)/me);
    % Decomposer ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    l = 1;    
    N = size_each_level(l);
    for n = 1 : N
        index_n = n + lb_level(l) - 1;
        index_n = sort_mass(index_n,2);
        vel(index_n,:) = rand(1,D).*(pos(sort_mass(ps,2),:)-pos(index_n,:))+pos(index_n,:);
%         vel(index_n,:) = rand(1,D).*(pub(index_n,:)-plb(index_n,:))+plb(index_n,:);
    end
    % Decomposer ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    % Consumer   ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    for l = 2 : num_FClevel-1
%         if((l == 2) & (mod(iter, 12) > 8))
% %         if(mod(iter, 12) > 8)
%             prob = rand();
%             dir_forage = (prob>=prob_invdir) - (prob<prob_invdir);
%         else
%             dir_forage = 1;
%         end
        N = size_each_level(l);
        level_mass = sort_mass(lb_level(l):ub_level(l),1);
        for n = 1 : N
            index_n = n + lb_level(l) - 1;
            index_n = sort_mass(index_n,2);
            if(rand < prob_forage)
%                 foraged_level = l + dir_forage;
                foraged_level = l + 1;
                index_prey = fix(rand*(ub_level(foraged_level)+1-lb_level(foraged_level)))+lb_level(foraged_level);
                index_prey = (index_prey>ub_level(foraged_level))*ub_level(foraged_level)+(index_prey<=ub_level(foraged_level))*index_prey;                
%                 index_prey = n + lb_level(foraged_level) - 1;
                index_prey = sort_mass(index_prey,2);
                % PSO
%                 vel(index_n,:) = pos(index_n,:) + (pos(index_prey,:)-pos(index_n,:)).*(2*rand(1,D)-1);
                % GWO
                vel(index_n,:) = ForageGWO(a, pos(index_prey,:), pos(index_n,:));
                % QGSA
%                 vel(index_n,:) = ForageQGSA(g, pos(index_prey,:), pos(index_n,:));
            else
                index_interact = RouletteWheel(level_mass,1);
                index_mate = index_interact;
                index_mate = index_mate+lb_level(l)-1;
                index_mate = sort_mass(index_mate,2);
                % PSO
%                 vel(index_n,:) = pos(index_n,:) + (pos(index_mate,:)-pos(index_n,:)).*(2*rand(1,D)-1);
                % GWO
                vel(index_n,:) = ForageGWO(a, pos(index_mate,:), pos(index_n,:));
                % QGSA
%                 vel(index_n,:) = ForageQGSA(g, pos(index_mate,:), pos(index_n,:));
            end
        end
    end
    % Consumer   ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    % Autotroph  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    l = num_FClevel;
    N = size_each_level(l);
    level_mass = sort_mass(lb_level(l):ub_level(l),1);
    index_best = sort_mass(ps, 2);
    for n = 1 : N-1
        index_interact = RouletteWheel(level_mass,1);
        index_n = n + lb_level(l) - 1;
        index_n = sort_mass(index_n,2);
        index_mate = index_interact;
        index_mate = index_mate+lb_level(l)-1;
        index_mate = sort_mass(index_mate,2);
%         % PSO
% %         vel(index_n,:) = pos(index_n,:) + (pos(index_mate,:)-pos(index_n,:)).*(2*rand(1,D)-1);
%         % GWO
%         vel(index_n,:) = ForageGWO(a, pos(index_mate,:), pos(index_n,:));
        % QGSA
       vel(index_n,:) = ForageQGSA(g, pos(index_mate,:), pos(index_n,:));
    end
    % Autotroph  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 
%     vel = ((vel>=vlb)&(vel<=vub)).*vel+(vel<vlb).*vlb+(vel>vub).*vub;
    pos = vel;
    pos = ((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*plb+(pos>pub).*pub;
    
    fitness = feval(fhd,pos',varargin{:});
    sort_fitness = [fitness' (1:ps)'];
    sort_fitness = sortrows(sort_fitness, 1);
    fitbest = sort_fitness(1,1);
    fitworst = sort_fitness(ps,1);
    gbestindex = sort_fitness(1,2);
    gbest = pos(gbestindex,:);
    gbestval = fitbest;
    
    if(bool_Plot)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
end

end
